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Abstract 

We derive new explicit expressions for the components of Moore-Penrose inverses of sym¬ 
metric difference matrices. These generalized inverses are applied in a new regularization 
approach for scattered data interpolation based on partial differential equations. The columns 
of the Moore-Penrose inverse then serve as elements of a dictionary that allow a sparse signal 
approximation. In order to find a set of suitable data points for signal representation we apply 
the orthogonal patching pursuit (OMP) method. 

Key words: Moore-Penrose inverse, Fourier transform, linear diffusion, orthogonal 
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1 Introduction 

Within the last years, there have been several attempts to derive new sparse representations 
for signals and images in the context of sparse data interpolation with partial differential 
equations. The essential key of these so called inpainting methods is to fix a suitable set of 
data points that admits a very good signal or image approximation, when the unspecified 
data are interpolated by means of a diffusion process. Finding optimal sparse data leads to a 
challenging optimization problem that is NP-hard, and different heuristic strategies have been 
proposed. Using nonlinear anisotropic diffusion processes, this approach has been introduced 
in and strongly improved since that time by more sophisticated optimization methods such 
as m- Inpainting based on linear differential operators such as the Laplacian is conceptually 
simpler, but can still provide sparse signal approximations with high quality; see e.g. Emms]. 

These recent results outperform earlier attempts in this direction that use specific features such 
as edges El [SI El or toppoints in scale-space m- 

In this paper, we focus on the one-dimensional discrete case. The goal is to provide new 
insights into this problem from a linear algebra perspective. In the first part we derive a formula 
for the Moore-Penrose inverse of difference matrices. This formula has already been proposed 
for finding generalized inverses of circulant matrices of size NxN with rank A—1 in [20] and for 
discrete Laplace operators considered e.g. in [HIITj. We generalize this formula to arbitrary 
symmetric NxN difference matrices of rank N —1 possessing the eigenvector 1 = (1,..., 1)^ 
to the single eigenvalue 0. This formula enables us to derive explicit expressions for the 
components of the Moore-Penrose inverse for special difference matrices that play a key role 
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in discrete diffusion inpainting. In particular, we study difference matrices that approximate 
the second and fourth order signal derivative, i.e., the one-dimensional Laplace operator and 
the biharmonic operator, with periodic resp. reflecting (homogeneous Neumann) boundary 
conditions. 

With the help of the generalized inverse of the difference matrices, we propose a new 
regularization approach that admits different optimization strategies for finding a suitable 
sparse set of data points for signal reconstruction. In particular, we employ the orthogonal 
matching pursuit (OMP) method as a conceptually simple and efficient greedy algorithm to 
construct the desired set of data points. 

The paper is organized as follows. In Section 2 we present the formula for the Moore-Penrose 
inverse of symmetric difference matrices and derive explicit expressions for the components of 
the Moore-Penrose inverse of difference matrices corresponding to the discrete Laplace and 
biharmonic operator using the theory of difference equations. Section 3 is devoted to a new 
regularization approach to the discrete inpainting problem. The columns of the explicitly given 
Moore-Penrose inverses of the difference matrices can be understood as discrete Green’s func¬ 
tions [8] and are used as a dictionary for the orthogonal matching pursuit algorithm (OMP). 
Since the original discrete inpainting problem (with the Laplace or biharmonic operator) is 
equivalent to a least squares spline approximation problem with free knots, our approach can 
also be seen as an alternative method to solve this nonlinear approximation problem. Finally 
we present some numerical results in Section 4. 


2 Explicit Moore-Penrose inverses for difference matrices 


Let L be an V X V matrix. We say that L is a difference matrix if it satisfies the following 
properties. 

(i) The matrix L is symmetric and of rank iV — 1; 

(ii) LI = 0 where !:=(!,..., 1)^ S K^, i.e., 1 is the eigenvector of L to the single eigenvalue 

0 . 

Matrices of this type typically occur when discrete approximations of symmetric differential 
operators with periodic or Neumann boundary conditions are considered. 

In particular, a discretization of the second derivative (one-dimensional Laplace) with pe¬ 
riodic boundary conditions yields a circulant matrix L of the form 


L = A = circ (-2,1,0,..., 0,1)^ := 


( -2 1 0 ... 1 \ 

1 -2 1 


0 


V 1 0 ... 1 -2 / 


while for Neumann conditions, we consider 


/ -1 1 0 

1 -2 1 

0 


L = B := 


Vo 0 


t,NxN 


( 2 . 1 ) 


0 \ 


-2 1 
1 -1 / 


cWxAf 


( 2 . 2 ) 


As a discretization of the fourth derivative (one-dimensional biharmonic operator) we obtain 
in the periodic case 

L = -A^ = circ (-6,4, -1, 0 ..., 0, -1, if 
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and in the case of reflecting boundary conditions 


L = -B^ 


All these matrices satisfy the properties (i) and (ii) given above. 

The following theorem generalizes the ideas of Eoiia for pseudo-inverses of circulant ma¬ 
trices of rank iV — 1, as well as of [min] for pseudo-inverses of the Laplacian. According to 
the usual definition, we say that L+ is the Moore-Penrose inverse of L if 

LL+L = L, L+LL+ = L+, 

and LL+ as well as L+L are symmetric. With these properties, the Moore-Penrose inverse is 
uniquely defined; see [Si- 

Theorem 2.1 For a difference matrix L G satisfying the conditions (i) and (ii) we 

obtain the Moore-Penrose inverse in the form 

L+= (L - rJ)-i + ;^J, (2.3) 

where r G M \ {0} can he chosen arbitrarily, and J is the N x N matrix where every element 
is set to 1. In particular, L"*" does not depend on the choice of t. Further, L’*' is symmetric. 

Proof. First, we show that L — rJ is a regular matrix for all r G K \ {0}. Since L 
is symmetric, there exists an orthonormal basis of eigenvectors vo,...,V 7 v-i for L, where 
vq = is the eigenvector of L to the single eigenvalue Aq = 0. Hence, Lvj = A^Vj, 

j = 1,... ,N — 1, with eigenvalues Aj G K \ {0} implies 

(L - rJ)vj = XjVj 

since l^Vj = 0 for all j = 1,..., — 1. Further, 

(L - rJ)l =-rJl =-riVl, (2.4) 


i.e., (L — rJ) possesses the N — 1 nonzero eigenvalues Ai,..., Xn-i of L and the additional 
eigenvalue —tN. Therefore, it is regular. In order to show that L+ in (2.3) is the Moore- 
Penrose inverse of L, we check the four properties that determine tie Moore-Penrose inverse 


uniquely. From (2.4) it follows that 


( I -—=-^1 


(2.5) 


and hence 

(L-tJ)-IJ 

Analogously, we simply derive J(L — tJ)“^ = 
JL = LJ = 0 we find 


tN 

-:^J from 1^(L 


rJ) = Using 


LL+L 


L(L-TJ)-i(L-TJ-krJ) 

L(I-k(L-TJ)-Vj) 



L, 


( 2 . 6 ) 
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and by = N3, 


L+LL+ = L+L (L - rJ)--^ + 


tN^ 

L+{L-t3 + t3) ( (L-rJ) 
1 


-1 


1 


= L+ 

= L+ 

= L+ 

= L+- 


r7V2 


(L-tJ)J + tJ(L-tJ)~i 


riV2 


I- 


— J2-—j+lj 

tN^ tN ^ N 


I-J 

N 


(2.7) 


N 


(L-rJ)- 


= L+-^ - 


N 


1 


+ 


1 


tN tN 


1 

j 


j J 


= L+. 


The equalities (2.6) and (2.7) particularly imply LL+ = L+L = (I — ;^J), i.e., L+L as well 
as LL+ are symmetric. Since L+ is uniquely defined, it does not depend on the choice of r. 
Finally, we have 


(L+)' =(lT)+=L+ 


which shows the symmetry of L+. 


□ 


We want to apply (2.3) in order to derive explicit expressions for the special difference 
matrices A and B in (2.1) and (2.2) as well as for — and —B2 that approximate differen¬ 
tial operators of second and fourth order, together with periodic or reflecting (homogeneous 
Neumann) boundary conditions. 


Theorem 2.2 The Moore-Penrose inverse of the circulant matrix A. in (2.1) is a symmetric 
matrix of the form 


( 


A+ =circ(aJ,a(^,...,a)()_J = 



a)v_i 

®)v-2 

a 

4 


a)v-i 


a-t 


<4 

a- 

+ 

_-i- 


+ 

N-l 

1 

to 




^ \ 

+ 

2 
+ 

3 




where 


< = 




j=0,...,A-l. 


Proof. By Theorem 


2.1 


with T = —1 we have A+ = (A -|- J)“^ — ]^J, where J is the 
matrix with every element one. Since A -|- J is circulant, also (A -|- J)“^ is circulant, and we 
assume that it is of the form (A -|- J)“^ = circ (6oj &i, • ■ •, ^ 7 V-i)- The proof of Theorem 
implies 


2.1 


A(A + J)-i = (A + J - J)(A + J)-i = I - J(A + J)-i = I - 


i.e., we obtain 


bj—i 2bj bjj — 1,...,A 1 


( 2 . 8 ) 
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and the boundary conditions b^-i — 26o + ~ ^Af -2 ~ 26jv-i + 6o = Further 

by (2.5) it follows that (A +J)“^l = ^1 and hence 


N-l 

E bj = 

3=0 


1 

N' 


(2.9) 


The inhomogenous difference equation (2.8) can be rewritten as 

0 


bj+i 


0 1 

-1 2 


6,-1 


' N 


With Yj+i := {bj,bj+i)'^, g := (0,-^)^, and 
order difference equation 

Vj+i = Gy 3 + s 

with the general solution 

y, =G^yo + ^G^-''g 


0 1 

-1 2 


j = l,...,7V-2. 


we obtain the linear, first 


where yo = (bo,bi)'^ needs to be determined by the boundary conditions and the additional 
condition (|2.9|). A simple induction argument shows that the powers of G can be written as 


G^ = 


-j + 1 j 
-j j + 1 


and we obtain 

bj = (1,0) yj 

= (1,0) 


- j + 1 j 

-j j + 1 


bo 

bi 


E(l^O) 


-j + 13+1 j -13 
V-j j-v+l 


' N 


= bo+j{bi-bo)- 


1 j{j - 1) 


N 2 

Finally, we determine bo and bi from 


N-l 


bN-i-2bo + bi = l- E = TV 

3=0 


( 2 . 10 ) 


( 2 . 11 ) 


while the second boundary condition is then automatically satisfied. We insert (2.10) into the 
first equation of (2.11) and obtain 

bo + {N- l)(6i - bo) - — - 26o + = 1 - 


i.e., 

“ 2N 

Hence, by (2.10), the second condition in (2.11) yields 

Af-l 

h. = Mh. ^ W - 

2N 


A - 1 1 1 

E E J - E = a 7 > 


3=0 


3=0 


3=0 


I.e., 




1 A2 - 1 


12 
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we find 


With bo = ^ 


N^-1 

12N 


, _ 1 I I jjN-j) 

^ iV2 12iV 2N ' 


We finally obtain from (2.3) (with r = — 1) 

1 




l-W" j{N-j) 

12N 2N 


Obviously A+ is symmetric and we have □ 

Further, we can directly derive the Moore-Penrose inverse for the discrete biharmonic op¬ 
erator —for periodic boundary conditions. 

Theorem 2.3 The Moore-Penrose inverse for the discrete biharmonic operator h = — A^ with 
A in (2.1) with periodic boundary conditions is a symmetric matrix of the form 

L+ = -(A2)+ =circ(c(J',c(^,...,c(!^_i) 

with 


c+ — 


(1 - A2)(A2-p 11) k{N - k){Nk--^2) 


-b 


fc = 0,...,A-l. 


720A 24A 

Proof. From (A^)+ = (A'*’)^ it follows that 

k N-l 

1=0 j=k+l 

Inserting ^ direct evaluation gives the desired result. Here one needs the 

Faulhaber formulas for sums of the form ft?-10. ^ where denotes the 

Bernoulli polynomial of degree f -b 1; see e.g. [T2]. □ 

As an interesting side information, we obtain the following trigonometric identities. 
Corollary 2.4 For N £ N we have 


^ COS(^) 

^ Sin2(t) 
cos(g#) 
sin^(^) 


- 1 


-2j(A-j), j=0,...,A-l, 


E 


(A2 - 1)(A2 + 11) 2j{N - j){Nj - f + 2) 


45 


j = 0,..., A - 1, 


y- sin(g^) 


and particularly 

N-l 


N-l 


V 3in(g^) 

h ’ 


j = 0,...,A-l, 


E 


- 1 


sin^(f) 


N-l 

E 


sin^^) 


(A2 - 1)(A2 -P 11) 
’ 45iv 


Proof. Since A in (|2.1|) resp. A^ are circulant matrices, they can be diagonalized by 

Af-l 


the Fourier matrix Fjv = ) 


diag with the characteristic polynomial a{z) ■= —2-\- z-\- ^ of A, i.e., 

a(^^) = _2 + 6-2-^=/^ + 62"^=/^ = _2 + 2 cos = -4 sin^ = 0,..., A - 1. 


l,fc=o 


, where wn := e 27i-*W, yye obtain A = FArAF)(r = 


6 


























Observe that a{uj^) are the eigenvalues of A resp. A. Considering the generalized Moore- 
Penrose inverse, we obtain the circulant matrix A+ = F^A+Fat with 


A+ := diag 0, 


-1 


-1 


-1 


4sin2(^) ’ 4sin2(^) ’ ’'' ’ 


Analogously, the Moore-Penrose inverse of A^ can be written as (A^)+ = F^(A+)^F 7 v. A 
comparison of the entries a'j (resp. Cj ) in the first column of of F^A+F^r, (resp. F^(A+)^FAr), 


N-l 


= -y 

N ^ 


-1 


1 me i ein 

,-ifc, ,feo _ ^ cos jy + z sm 


k=l 


N . 


4sin2(^) 


N ^ 




4sin^(f) 


resp. 




A^—1 2'Kkj I • • 2'Kkj 

1 COS —+ I Sin-^ 


= -F-- 

N ^ 


N 


N 


16sin4(^) 

with the formulas found for A+ (resp. (A^)+) in Theorems 2.2 and 2.3 yields the assertions. 

□ 


Let us now consider the case of reflecting boundary conditions and derive the Moore-Penrose 
inverse for the matrix B in (2.2) and for —B^. 


Theorem 2.5 The Moore-Penrose inverse of the matrix B in (2.2) occurring by discretization 
of the Laplacian in the case of homogeneous Neumann boundary conditions is a symmetric 
matrix of the form B+ = {b'^^.)^ifJ:Q with 

(A-l)(2iV-l) , , j(j + l) Hk + l) 

=-6iV-^ - 

for j > k. 

Proof. 1. Let now A = A 2 N = circ (—2,1,..., 1)^ S I 


2N 


,2Nx2N 


be the circulant matrix as 


in (2.1) but of order 2N. We consider the four N x N blocks of A 2 N, 

A2N = 


Ai A 


Af 

0 


where Aq and Ai are symmetric Toeplitz matrices of size N x N. In particular, Ai = Af has 
only two nonzero entries. Further, we denote by A^r = circ (a^", a^,..., the Moore- 

Penrose inverse of A 2 JV, where ~ 4A^) -|- ^k(2N — k) by Theorem 

N X N blocks of A^j^ be denoted by 


2.2 


Let the 


A+ — 

^2N — 


K 


At^ 

A + 

“^0 


Observe that also A,J" and A((~ are Toeplitz matrices. Further, let I denote the counter identity 
of size N X N, i.e., 

f M 


1 


V 1 




Then we simply observe that the matrix B in (2.2) can be written as 

B = Aq -|- lAi. 
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We show now that B+ := Aj + lA^ is the Moore-Penrose inverse of B, where Aj and A^ 
are the partial matrices of A^^ as given above. We know that 


^2N-^2N-^2N — 




Aq 

Ai 



AT 

Aq 


which leads to 


AoA+Ao + AfA+Ao+Ao(A+)^Ai+AfA+Ai = Aq, 

AiA(j“Ao + AqAj'^Ao + Ai(A^)'^Ai + AoA(j“Ai = Ai. 

Using these equalities, and taking into account the symmetry of Aq, Ai and A,^ as well as 
that ITI = for any Toeplitz matrix T, a short computation gives 


BB+B = (Ao+iAi)(A++iA+)(Ao + iAi) 

= AoA+Ao + AoA+iAi + AqIA+Ao + AoIA+iAi + iAi A+Aq + iAi A+IAi 

+iAiiA+Ao + lAiiA+IAi 

= AoA+Ao + iA^(A+f Af + iA;^A+Ao + Ao(A+f Ai + IAiA+Aq 
+AT (A+fAi + AfA+Ao + iAi(A+)^Ai 
= (Ao+iAi)=B 


and analogously B+BB+ = B+. Further, using the symmetry of A 2 AfA^^ we observe that 
(A„A+ + Af A+) is symmetric yielding 

BB+ = (Ao + iAi)(A++iA+) 

= AoA+ + iAiiA+ + AoiA+ + iAiA+ 

= ( Aq A(j" + AT AT ) + i( A(^ A^ + Ai A|^) 

= (AoA+ + AfA+f+ (AjA++AiA+fi 
= (A+fAT + (AtfA, + (ATfAoi + {A+fATi 
= (A++iA+)^(Ao + iAif = (BB+)^ 


and analogously B+B = (B+B)^ 


2. Using the formulas for the (2N x 2N) matrix AJv from Theorem 

B+ — (h~^ \N —1 A+ — (fj+ \ 2 N —1 — / + \2Af—1 

“ \°i.kU.k=0 ^2N — \^j,k)j,k=Q “ l®0-fc)mod 2Af U.fc=e 


2.2 


we now obtain 


0 


Kk 



+ d 


+ 

2N-l-j,k 


^{j-k)mod 2N 


+ CL 


+ 

{2N—j — k — l)mod 2N‘ 


For j > k this gives 


U = 2(1-(2A)^) (j - k){2N - j + k) i2N-j-k-l)ij + k + l) 

'3,k - 24N ■■ 4iV 4A 

(A-1)(2A-1) . jU + 1) kjk + l) 

6N 2N 2N ' 


□ 


Theo rem 2.6 

B in i2.2\ is a 


The Moore-Penrose inverse for the discrete biharmonic operator h = —B^ with 
symmetric matrix of the form L+ = —(B^)+ = {d'^jf)fk=o 


(A" - 1)(4A^ - 1) [j{j + l) + k{k + l)f j{j + l)k{k + l) 

180A 24A 6N 

+ j iUU + 1) + Hk + 1)] - [j{j + 1) + nk + 1)] 


for j > k. 





















Proof. From Theorem 
j < k, where 


2.5 


Pj,k ■ — 


it follows that = j3j^k + j for j > k and b^f, = Pj^k + k for 
(fV-l)(2iV-l) j{j + l) k{k + l) 


■ 6N 2N 

Observing that (B^)+ = (B+)^, we obtain for j > k 

k j 


2N 


N-l 


+ j)(/3£,fc + fc) + ^ iPj,e + j){l3i,k + ^) + ^ {Pj,e + ^){Pe,k + ^)- 

e=o i=k+i ^=i+i 


Inserting the value Pj/, some algebraic evaluations give the desired result. 


□ 


3 Sparse signal approximation by regularized diffusion 
inpainting 

Using the above explicit representations of the Moore-Penrose inverses of difference matrices, 
we want to derive a new algorithm for sparse signal approximation based on a regularization 
of a discrete diffusion inpainting method. 


3.1 Regularized diffusion inpainting 

We consider now the following approach for sparse signal representation. For a finite signal 
f = we want to find a good approximation u = {uj)fSQ G using only 

a small amount of values fk, with k G F C := {0,...,iV — I}. Here, F denotes an index 
set with a small number |r| of entries, i.e. |F| N. In order to construct u, we apply an 
inpainting procedure. We want to approximate f by finding a solution of 


(Lu)fc = 0 forfce {0, ...,fv-i}\r, 
Uk = Ik for k GT, 


(3.12) 

(3.13) 


where L denotes a symmetric difference matrix satisfying the conditions (i) and (ii). Here, 
L can be seen as a smoothing operator, as e.g. the discretization of the second derivative 
(Laplacian) or the fourth derivative (biharmonic operator) with periodic or reflecting boundary 
conditions. 

Let the vector c := (cq, ..., cn-i)'^ be given by 


r 1 iGT 
\ 0 i e H\F, 


Further, let C := diag (c) be the diagonal matrix determined by c, and let In denote the 
identity matrix of order N. 

No we can reformulate the problem as follows. Find an index set F (or equivalently an index 
vector c) with a fixed number of values |F| ^ N, as well as a sparse vector f = (/o,..., Jn-i)^ 
with /fc = 0 for k G ft\r, such that u can be recovered from f as 


C(u- f) - (Itv - C)Lu = 0 


(3.14) 


and 11 u — f ||2 is (in some sense) minimal. 

In the two-dimensional case with homogeneous Neumann boundary conditions, Mainberger 
et al. m considered a specific greedy algorithm: To determine a suitable set F, they used a 
probabilistic sparsification followed by a nonlocal pixel exchange. Afterwards, the tonal (i.e. 
greyvalue) data f{k) for k gT are optimized with a least squares approach. Here we propose a 
different strategy that simultaneously optimizes spatial and tonal data. Equation (3.14) yields 


Lu = C(u-f-hLu). 
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Since C is a sparse diagonal matrix with only |r| nonzero diagonal entries, it follows that 
g := C(u — f + Lu) is a sparse vector with |r| nonzero entries. Therefore, we reformulate our 
problem: We aim to solve 

Lu = g, 

where the sparse vector g needs to be chosen in a way such that the error ||u — f ||2 is small. 
Since L is not invertible we apply the Moore-Penrose inverse of L to obtain 


u = L+g. (3.15) 

From l^L = 0^, we find that l^Lu = l^g = 0. This solvability condition is satisfied if 
l^C(u — f — Lu) = 0 is true. Further, we also have 1^L+ = 0^, i.e., each solution u = L+g 
in ( 3.15| ) has the mean value 0. Since we are interested in a sparse approximation u of f, we 
therefore also assume that also f has mean value 0. In practice, we just replace f by f° = f — /I 
with 

1 1 

3=0 

and store / elsewhere. 


Remark. The regularization approach described above to solve the inpainting problem can 
also be interpreted as a solution using the concept of discrete Green’s functions; see [8] for the 
two-dimensional case. 


3.2 Orthogonal matching pursnit approach 

With the regularization approach ( |3.15 ), the original task of finding a suitable index set F 
and the corresponding optimal values /fe, fc S F can be rewritten as the following optimization 
problem: Find g S 


r,N 


such that 

l|L+g- f ||2 min subject to ||g||o = |r|. 


where ||g|io denotes the number of nonzero coefficients in the vector g. Let us denote the 
columns of the matrix L'*' by ag,..., ajv-i. We now rewrite the optimization problem in the 
form, 

||zz|| 2 ^-min subject to L+g-f zz = f, ||g||o = |r|. 

Equivalently, with F denoting the index set of the nonzero entries of the sparse vector g and 
|r| = ||g|lo "C N, we have f in the form 


f = gfc afc + 

where f is approximated by a sparse linear combination of the columns a^., k € F. We apply 
now the orthogonal matching pursuit as a greedy algorithm for finding a suitable subset F as 
well as the sparse vector g that determines u. Here, {ag,..., ajv-i} is the dictionary for our 
OMP method in the Hilbert space This algorithm is known to work very efficiently also 
when the original signal f is not exactly sparse in the considered dictionary, see e.g. [2T| . 

The OMP works as follows. In a first step, we determine the index ki G {0,..., — 1} 

such that the column a/jj^ correlates most strongly with f, i.e. 

, l(f,ak)| 

ki = arg max -r, 

k=o,...,N-i (ak,ak) 

where (f, a*,) = f^aj, is the standard scalar product of the two vectors f and a^. 

In the next step, we determine the coefficient gk^ such that the Euclidean norm ||f — gk^ ■ aki II 2 
is minimal, i.e. gk^ = (f, afcj^)/||ai,J||, where ||ai;J |2 denotes the Euclidean norm of a^^. 

Now we consider the residuum ri = f — gki&ki and proceed with the next iteration step, 
where f is replaced by ri. For all further iteration steps we add an update of the coefficients. 
The algorithm can be summarized as follows. 
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Algorithm 3.1 (OMP) 

Input: Dictionary {ap,..., a 7 v-i} o/K^, rp = f S (with l^f = 0^, g = 0, |r| <C N. 
Forj = l,..., |r| do 

1. Determine an optimal index kj such that a^^ correlates most strongly with the residuum 


rj_i, i.e. 


kj = arg max 




fc=o,...,JV-i (afe,afe) 

2. Update the coefficients , • ■ •, gUj such that ||f — 9ki ^ki II 2 is minimal and set Vj = 
end (do) 

Output: gi' = (gfe.)lL'i 


As proposed in the algorithm, for a given number |r| ^ N, we may just take |r| iteration 
steps. For improving the method, one may also take some more steps, and afterwards decide 
by a shrinkage procedure, which indices are kept to determine F. 

Remarks. 1. For numerical purposes, we normalize the columns ap,... ,ajv_i in a prepro¬ 
cessing step before starting the OMP algorithm. Using L+ based on the Laplace operator 
or the biharmonic operator with periodic boundary conditions, see Theorems 2.2 and 2.3 all 
columns have the same Euclidean norm, and we obtain for example 


N-l 

E( 


at? = 


N-l 

E 

fe =0 


1 - 
12N 


k{N-k) 

T/V 


{N^ + 11)(A2 - 1) 

720A 


for the Moore-Penrose inverse in Theorem |2.2| Note that the columns occurring in the Moore- 
Penrose inverses in Theorems |2.5| and |2.6| for reflecting boundary conditions do not have the 
same Euclidean norm, such that the normalization in the OMP algorithm is crucial. 

2. Observe that the considered regularization does not longer exactly solve the original 


inpainting problem (3.12)-(3.13) resp. (3.141, since we cannot ensure the solvability condition 
l^g = 0. One way to enforce this condition is by incorporating it as additional constraint, see 
Section 4. 


3.3 Relation to spline interpolation with arbitrary knots 

Finally, we want to point out the close relationship between the proposed inpainting algorithm 
considered above and the nonlinear spline approximation problem with arbitrary knots in order 
to find a sparse approximation of f. For that purpose, we consider the inpainting problem 


(3.14) for the special case L = B in (2.2) corresponding to the Laplace operator with reflecting 


(homogeneous Neumann) boundary conditions. 


Cu- (Iat - C)Bu = Cf. 


(3.16) 


Assume that the index set of 1-entries of the index vector c is F = {£ 1 ,... ,£^1, where 0 < 
(-1 < (-2 < ■ ■ ■ < < N — In this case, the matrix M := C — (Ijv — C)B is an invertible 

block diagonal matrix with L -|- 1 blocks Mj, where the inner blocks are of the form 


M,= 


/ 1 0 

1 -2 

0 

VO ... 


1 -2 

0 0 


0 \ 
0 


1 
1 / 


i+i-L+i)x(L+i-L+i) 


£ = 1,...,L-1, 
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and the boundary blocks 


Mn = 


1 

-2 1 

0 c 


1 -2 

1 

0 c 




II 




0 c 

1 

-2 1 


0 c 

1 

-2 1 


are of order + \ resp. N — II- Observe that the boundary blocks can degenerate to Mg = 1 

^ ^ ' for = 1, and analogously, = 1 for = TV — 1 and 


for £i = 0 and Mg = 


0 1 


Mr, = 


1 0 
1 -1 


for = N — 2. The inner blocks degenerate to = I2 for ^i+i - + 1- 


Now the system Mu = Cf in (3.16) can be solved by solving the partial systems 


where 


'-^3 

Cl 

Cl 

yd) 

f(i) 


diag(l,0,...,0,l) e k(^.+i-^j+i)x(^,+i-^,+i) £ = !,.. 

diag(0,...,0,l) e K(^i+i)x(fi+i)^ 
diag(l,0,...,0) G k(^-^i^)x(JV-^i.)^ 
j = 0,...,L 

( 4 , 0 ,..., 0 , 4 ^,f j= 0 ,...,L. 


For the inner blocks we observe that = (4., 0,..., 0, )^, and it can be simply 

verified that 


UE- + k = 


{rij - k)fi + kfi 


rij := lj+\ — iji k = 0,... ,nj. 


For the boundary blocks, the constant vectors 

u(°)=4 4, 

solve the partial systems. As a global solution, we hence obtain a discrete linear spline ap¬ 
proximation u of f with knots given by the index set F resp. c. Thus, the inpainting problem 
considered in (3.16) is equivalent to the nonlinear approximation problem to find optimal knots 
ij (out of the set {0 ,..., A^— 1 }) and corresponding values 4 > such that the obtained discrete 
piecewise linear spline vector 


u = ((u(°))^,...,(u(^))^)^ G 


cAf 


minimizes the norm ||f — u|| 2 . 

As we will see in the next section, our regularization approach using the OMP method 
can be thus also seen as a method to find an solution for the discrete spline approximation 
problem. 

Remarks. 

1. Analogously, it can be shown that the inpainting algorithm with the operator —corre¬ 
sponding to the biharmonic operator with reflecting boundary conditions is related to nonlinear 
discrete spline approximation with cubic splines with arbitrary knots. 

2. Sparse data approximation using splines with free knots has been already extensively 
studied in the literature, see e.g. [IIlIlsllISllI] and references therein. Splines with free knots 
have been shown to perform significantly better than splines with fixed knots. For a discussion 
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in the context of signal approximation we refer to [7]. The question to find the optimal knot 
sequence leads to a global optimization problem that is nonlinear and non-convex. Therefore, 
besides expensive global optimization strategies, there exist also heuristic strategies, where one 
starts either with only a few knots and iteratively adds further knots to improve the quality 
of the approximating spline, or conversely, starts with a large number of knots and iteratively 
removes knots that are less important for the approximation error. Our approach belongs to 
the first group. 


4 Numerical examples 

We apply now the regularization of the inpainting approach proposed in Section 3 to sparse ap¬ 
proximation of vectors. In particular, we want to use the Moore-Penrose inverses of the discrete 
Laplacian and the discrete biharmonic operator with either periodic or reflecting boundary con¬ 
ditions that have been explicitly computed in Section 2. Our previous considerations imply 
the following simple algorithm. 


Algorithm 4.1 (Sparse approximation) 

Input: {ao,..., aAr_i} columns of the matrix L+, 

f S |r| -|- 1 number of values used for sparse approximation 


1. Compute f := ^l^f and := f — /I. 


3.1 


2. Apply the OMP-algorithm 
o/fO. 

Output: f, g'" = determining a sparse approximation oft of the form 


Tl 


to with |r| steps to obtain an approximation ^ g^- 

Z = 1 


|r| 

u = + afe,. 

i=l 


We have applied this algorithm using the columns of the Moore-Penrose inverse of L for 
L e {A, B, -A2, -B^} as computed on Section 2. Starting with a vector f of length N = 256, 
we spend |r| = 13 values for the approximation of f, i.e., about 5 % of the original vector size. 
In Figure 1, we present the approximation results together with their error. As is common in 
signal processing, we measure the error in terms of its peak signal-to-noise ratio (PSNR). For 
a maximnm possible signal value of M (in our experiments: M = 1), the PSNR between two 
signals u = and f = (/j)(l“Ms given by 

PSNR(u,f) := 10 1og„(jj^) (4.17) 

where MSE denotes the mean squared error 


N-l 

MSE(u,f) = -^(wg-/g)2. (4.18) 

i=0 

Note that higher PSNR values correspond to better approximations. Since the signal vector f 
in Fig. 1 is not “smooth”, the inpainting approximation with the discrete Laplacian provides 
better results than the discrete biharmonic operator. 

In the procedure described above, we have not implemented the “solvability condition” 

l^g = 0 (4.19) 
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homogeneous Neumann periodic boundary 

boundary conditions conditions 


Laplace operator 


Biharmonic operator 




PSNR: 25.78 


PSNR: 24.39 




PSNR: 25.58 


PSNR: 25.30 


Figure 1. Sparse approximation of a vector f € (red line) by a vector u (blue line) that is a linear combination 
of 13 columns of the pseudo-inverse L+ according to Algorithm 4.1. For L+, two different operators and two 
boundary conditions are considered. 


that we had derived from Lu = g and l^L = 0^ in Section 3. Therefore, the obtained approx¬ 
imation results in Figure 1 (top) for the Laplacian case appear to have piecewise quadratic 
form, since the columns of and B+ in Theorems 2.2 and 2.5 are quadratic in j. In order 
to solve the original inpainting problem, we need also to incorporate the constraint l^g = 0. 
This can be done by employing the additional equation 9ki = 0 in the last iteration step 
of the OMP algorithm when updating the coefficients gt^,. ■ ■, by solving the least squares 
minimization problem. 


The results of the constrained sparse approximation, i.e., employing also the solvability 
condition, for the four cases L € {A, B, —A^, —B^} are presented in Figure 2. Here we have 
used the same vector f S and |r| = 13 as in the first experiment. In this case, we observe 
that the inpainting approach with the Laplacian provides a piecewise linear approximation. 


Both figures illustrate that the choice of the boundary conditions can have a fairly strong 
impact on the approximation result, even in locations that are far away from the boundaries. 
In 3 out of 4 cases, homogeneous Neumann boundary conditions gave better approximations 
than periodic ones. This is to be expected for a non-periodic real-world signal. 
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Figure 2. Sparse approximation of a vector f € 


(red line) by a vector u (blue line) being a solution of the 


inpainting problem with different operators and boundary conditions. In contrast to Fig. 1, the solvability condition 
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